Exact rate calculations by trajectory parallelization and twisting 
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A sampling procedure to compute exactly the rate of activated processes arising in systems at 
equilibrium or nonequilibrium steady state is presented. The procedure is a generalization of the 
method in [A. Warmflash, P. Bhimalapuram, and A. R. Dinner, J. Chem. Phys. 127, 154112 (2007); 
A. Dickson, A. Warmflash, and A. R. Dinner, J. Chem. Phys. 130, 074104 (2009)] in which one 
performs simulations restricted into cells by using a reinjection rule at the boundaries of the cells 
which is consistent with the exact probability fluxes through these boundaries. Our generalization 
uses results from transition path theory which indicate how to twist the dynamics to calculate 
reaction rates. 
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Introduction. The main objective of this work is to 
revisit and extend in scope the nonequilibrium sampling 
procedure to compute steady state probability distribu- 
tions proposed in Rcfs. 1,2'. Specifically, we show how this 
procedure can be modified to calculate exactly certain 
dynamical quantities such as the rate of reactions occur- 
ring in arbitrary equilibrium or nonequilibrium systems 
at statistical steady state. This is done by exploiting re- 
sults of transition path theory (TPT}2<i^i^ which indicate 
how to twist the dynamics of a given system to calculate 
the reaction rate between a given reactant and product 
state. The method proposed in this note can be seen 
as a generalization to arbitrary nonequilibrium processes 
at statistical steady state of the milestoning procedure 
with Voronoi tessellation proposed in Ref. [gI by building 
upon the original works in Refs. This generaliza- 

tion makes the procedure more expensive computation- 
ally, but it permits to relax completely the assumptions 
made in milestoning - these assumptions were discussed 
in detail in Ref. Jii. Our method can also be viewed 
as a generalization of the transition interface sampling 
(TIS>iiii^ii^ and forward flux sampling (FFS)^^^^ meth- 
ods in which arbitrary sets of interfaces can be used that 
do not have to be placed in monotone succession. 

The remainder of this note is organized as follows. 
First we revisit from an original perspective the nonequi- 
librium sampling procedure of Refs. 1,2. Next, we show 
how this procedure can be modified to calculate reaction 
rates exactly using TFT. We then compare our proce- 
dure with the Markovian milestoning method proposed 
in Ref. 'd and with TlSiiii^ii^ and FFgi^'-i^. Finally we 
illustrate our procedure on a simple example. In terms 
of notations and assumptions, we will denote by z the 
location of the system in its state-space f2 C M'' (e.g. it 
could be the positions and velocities of all the atoms in a 
molecular system, in which case z ~ {x, v) and d ~ 6n ii 
n is the number of atoms). The specifics of the dynamics 
of the system are not important except that we assume 
that (i) its evolution is Markovian and (ii) it is ergodic 
with respect to a probability density function which we 



denote by g{z). Notice that we do not require detailed 
balance, i.e. g{z) is associated with a nonequilibrium sta- 
tistical steady state in general. 

Restricted sampling with flux matching. The methods 
m Refs.[l|ii are based on a factorization of the dynamics 
in which one artificially constrains the system to evolve 
in a set of cells partitioning state-space in a way that 
(i) does not bias the dynamics inside the cells and (ii) 
is consistent with the exact probability fiuxes in and out 
of these cells. In other words, the procedure guarantees 
that a true unconstrained trajectory of the system can 
be reconstructed exactly by patching together in some 
appropriate way the pieces computed in the cells. These 
pieces can be calculated in parallel in each cell, hence the 
name trajectory parallelization. The method in Ref. is 
restricted to equilibrium systems and exploits the time- 
reversibility of the dynamics. The method in Refs. [1113 
is more costly but it works also for nonequilibrium sys- 
tems. To explain how the latter method works, we first 
recall how to construct the cells using the Voronoi tes- 
sellation associated with a given set of generating points 
or center a^i^^ . Denoting these centers hy Za € ^ C M.'^, 
with a = 1, . . . , A, the Voronoi cell Ba associated to Za 
contains all the points that are closer to Za than to any 
other center, i.e. (see Fig. [1] for an illustration) 



{z e fl : ||2;-2;a|| < ||2;-2;/3|| for aU /? 7^ a}, (1) 



where || • || is some appropriate norm (e.g. the Euclidean 
norm in which case WzW^ — ^f)- 

If we were to generate an infinitely long trajectory of 
the system, z{t) with t > 0, this trajectory would keep 
going in and out of the cells Ba by crossing the edges be- 
tween these cells. Out of this trajectory we could there- 
fore generate an ensemble of exit-entry points, i.e. those 
points on the edges of the cells at which the trajectory 
goes from a cell Ba into a neighboring cell Bp. By the 
Markovian assumption, these points are all we need to 
generate an exact sample of trajectories inside each of 
the cells Ba simply by starting trajectories at the points 
leading into that cell and running these trajectories for- 
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FIG. 1: Representation of a portion of the Voronoi tessellation 
in two- dimension, where point Za is the generator of cell Ba 
and point the generator of cell Bj3. Pieces of restricted 
trajectory inside Ba are shown in black and inside Bp in grey. 
Each time the trajectory in one cell attempts to exit the cell, 
the exit point position (depicted as an open circle) is stored, 
and used to reinitialize the trajectory into the neighboring 
ceU. 



ward in time until they exit the cell. The procedure in 
Refs. [HIl is a way to generate these pieces of trajecto- 
ries inside the cells without having to compute the en- 
try points beforehand from a long unbiased trajectory 
but rather by generating them on-the-fly. To understand 
how this is done, imagine that we associate an indepen- 
dent copy, or replica, of the system to each cell Ba- Let 
us denote the instantaneous position of these replicas by 
G Ba-, a = 1,...,A. Even if we start Za{t) in- 
side Ba, sooner or later this trajectory will try to exit 
Be and go to another cell. When this happens, we store 
the exit point on the boundary, put the trajectory on 
hold and wait until a trajectory in one of the neighbor- 
ing cells makes an attempt to exit this cell by crossing 
the boundary with Ba- We then take this crossing point 
and use it to reinitialize the trajectory in B^. This is il- 
lustrated in Fig. [1] At the beginning of the simulation we 
do not have enough re-entry points, so many of the repli- 
cas may be on hold. But as the simulation goes on, we 
can build databanks of re-entry points from any Bp into 
any Ba (assuming that these two cells have a common 
boundary - otherwise the databank is trivially empty), 
which contain the last X points by which the trajectory 
tried to escape from cell Bp and enter cell B^ ■ This way, 
each time a trajectory tries to exit a cell Ba, we can im- 
mediately pick a re-entry point into that cell from the 
appropriate databank and continue the trajectory from 
that point without having to put it on hold. 

The only remaining issue we need to take care of in 
order to make complete the procedure outlined above is 
how to pick the re-entry point into Ba among the data- 
banks on the various edges leading into Ba- By con- 
struction, the re-entry points in the databanks on each 
edge are unbiased samples of re-entry points conditional 



on the trajectory entering by that edge. But there are 
several edges by which the trajectory can enter a given 
cell, and to introduce no bias we need to pick the edges 
with the proper probability of re-entrance by that edge. 
To see how this can be done, imagine that, as we run the 
simulation in each cell Ba, we compute an estimate of 
the effective rate of exit out that cell and into Bp via 



Na,p 
Ta 



(2) 



where Na^p is the total number of times the trajectory 
hit the boundary between Ba and Bp (which is also the 
number of times the trajectory in cell Ba had to be re- 
injected into that cell from a re-entry point) and Ta is the 
total simulation time in cell Ba (i.e. the total time the 
trajectory in cell Ba was running). If we then denote by 
tTq, the probability to find the unbiased trajectory inside 
cell Ba at statistical steady state, i.e. 



Q{z)dz, 



we see that tTq, and Va p are related as 



A 

E 

13=1 
/3/« 



T^pVp^a 



A 

E 

(3=1 
f3^a 



T^a'^a,f3, 



Y^TTa^l. 



(3) 



(4) 



The first equation in ^ simply expresses that, at sta- 
tistical steady state, the total probability flux into Ba 
(which is the term at the left hand-side of the first equa- 
tion in (HI)) must be equal to the total flux out of Ba 
(which is the term at the right hand-side of the first 
equation in (|3])). The second equation in ^ is simply 
a normalization condition for the probability which fol- 
lows from Ea^a = T^als^ g{z)dz = J^g{z)dz = 1. 
We stress that Q does not require detailed balance 
(i.e. TTpvp^a 7^ T^aVa,i3, possibly) and so it holds even for 
nonequilibrium processes provided only that they are at 
statistical steady state. Having calculated Va.p from ^ 
and TTq, from ([4]), we then have an estimate for the prob- 
ability flux from Bp into Ba- T^pvp^a- Consistently, the 
probability that the trajectory enters cell Ba by coming 
from Bp is simply: 



dB^ndB^ 



T^pVp^a 



(5) 



if Ba and Bp have a common edge and Pas^nOBc = 
otherwise. This expression gives us the desired probabil- 
ity to pick the edge dBp n dBa for re-entry into Ba - 

Summarizing, the algorithm to perform simulations re- 
stricted inside the cells for systems at nonequilibrium 
steady state is as follows: 

(1) Denoting by Za{t) E Ba the current state of the 
repfica in cell Ba, let 2;* be the state of the sys- 
tem produced from Za{t) after one timestep At by 
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a standard (i.e. unrestricted) integrator for the sys- 
tem (e.g. velocity- Verlet if the system is a molecular 
dynamics (MD) system). If 2;* G 5^, set 

z^{t + At) = zl. (6) 

Otherwise, if 2;* G Bp with (3 ^ a: 

(i) Store the point 2;* in a databank of entry 
points from Ba into Bp. In other words, as- 
suming that the databank contains already m 
points 2;^ p with fc = 1,2, ... ,m, set 2;™^^ = 

(ii) Update ^a,/3 via ([2|), tTq via (|4|) and Pas^jnas^ 
via (0). 

(iii) Select an edge dB^f n of B^ with proba- 
bility proportional to VgB^j,ndBc- 

(iv) Pick a point 2;^, with uniform probability 
from the databank of re-entry points on edge 
dBfj' n dBa , and set 

z^it + At) = z%^„ (7) 

(2) Go to step (1) and iterate to collect statistics. 

In essence, the algorithm above is the same as the one 
proposed in Ref. [HU, though the two differ in the details. 
For instance, we use a single set of cells instead of two (in 
Refs. 1,2 two staggered sets were used to ensure stabil- 
ity, but we observed no such stability problems with the 
procedure above). Besides a set of trajectories, one of 
the output of the procedure is to give the probability tTq 
to find the system in cell Ba at statistical steady state 
(see Indeed the computation of tTq, was the main 

objective in Refs.[l|H. Below we will show how to extract 
more from the procedure, in particular reaction rate in- 
formation, by appropriate modifications. Before getting 
there, however, let us make a few remarks about the al- 
gorithm above. In step (iv) we have assumed that the 
list of entry points from Bp' into Ba is not empty and, 
as already mentioned above, this may not be true at the 
beginning of the simulation (there might have been no 
collision with that edge up to that time). In that case 
we have to put the simulations in some of cells on hold 
at the beginning until we get proper re-entry points and 
start building databanks on the edges. Note that, at any 
given time, these databanks could each contain the last 
re-entry point only, though in that case we may run out 
of points again and have to put replicas on hold. Thus 
it is safer to always keep as many points in the data- 
banks as memory allows. However, the databanks do 
not have to be enlarged indefinitely - at worst we trade 
memory for CPUs since, the smaller the databanks, the 
higher the probability that one replica will have to be 
put on hold for some time. Also note that the quantities 
used to evaluate the probability in Eq. ([5]) have to be 
computed on-the-fly, and need information from all the 
cells. This, again, may lead to problems at the begin- 
ning of the simulations, when the statistics for Na.p is 



not accurate enough. To overcome this problem we set 
TTa — est Va = 1, . . . ,N at the beginning when statis- 
tics is insufficient to solve Eq. Using more educated 
guesses is possible too. Also, it is worth noting that i^ct,/3 
and TTa can be monitored on-the-fly to assess their con- 
vergence as a function of the length of the simulation, 
and the actual simulations of the replicas in each cell can 
be performed in parallel; only the re-entry events require 
communication. Finally, we should stress that the pro- 
cedure above relies on the ability to count the successive 
points at which a trajectory crosses the edges of the cells. 
This may lead to difficulties if the dynamics is governed 
by a stochastic differential equation, in which case these 
crossing points may form a fractal set. It leads to no 
difficulty, however, if some components of the trajectory 
are smooth and the cells are defined accordingly. We will 
come back to this issue later in the illustrative example 
section. 

Rate ealculation. Let us now come to the question of 
how to compute the reaction rate between a reactant and 
product state, which we identify as two disjoint sets in 
the system's state-space denoted as ^4 C O and B C 
respectively. As mentioned earlier, this calculation will 
be done by twisting the dynamics in some appropriate 
way that is dictated by the results of TPT. We begin by 
introducing the relevant objects that we will consider. If 
we assume again that we have at our disposal an infinitely 
long trajectory, z{t) with t > 0, this trajectory will go 
back and forth between A and B as time goes on, and 
we can split this trajectory into two pieces, depending 
on whether it visited last A or B. This construction is 
illustrated in Fig. [21 where the trajectory is shown in red 
if it visited A last and in black if it visited B last: if the 
trajectory visited A last at time t, we will say that it is 
assigned to A at time t, if it visited B last at time t, we 
will say that it is assigned to B at time t. 




FIG. 2: Schematic representation of a piece of ergodic tra- 
jectory visiting the two sets A and B. The pieces of this 
trajectory for which the last visited set was A are depicted 
in red, and those for which the last visited set was B are 
depicted in black. 

Based on this assignment, we can introduce the fol- 
lowing quantities. If iVj ^ denotes the total number of 
times the trajectory went from being assigned to A to be- 
ing assigned to B during the time interval [0, T] (i.e. the 
number of times it switched from red to black in Fig. [2]), 
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we set 



lim 



T ' 



(8) 



This quantity gives the average frequency at which the 
trajectory goes from Aio B which, because the system is 
at steady state, is also the same as the average frequency 
at which the trajectory goes from B io A since iVj^ = 
iV^^ asymptotically. The rate in ^ is referred to as 
the rate of the reactive trajectories in TPT, hence the 
subscript R. The rate of reactive trajectories vn should 
not be confused with the two reaction rates from Aio B 
and B io A defined respectively as 



kA,B = lim 

T— cx 



Ta ' 



^B,A 



lim 



NIb 
Tb ' 



(9) 



where Ta and Tb are, respectively, the total times dur- 
ing which the trajectory was assigned to A or i? in the 
interval [0, T] (i.e. the total times the trajectory is red 
or black in Fig. [2]). If A and B are metastable (i.e. if the 
trajectory commits to each of these sets and looses mem- 
ory of its past before going back to the other set), kA,B 
and kB,A are the rates that enter the phenomenological 
mass-action law describing how the populations in A and 
B evolve in time. Notice that is related to kA,B and 
kB,A as 



kA, 



VR_ 

pa'' 



^B.A 



(10) 



where pA is the fraction of time the trajectory is assigned 
to A and pB the fraction of time it is assigned to B: 



PA 



lim — ^ , 



PB 



lim — 



(11) 



By definition pa < 1, < 1 and pa + Pb = 1- The 
quantities i^n, kA,B, kB,A, PA and pB are the ones we 
will show how to compute exactly. To see how this can 
be done, it is useful to give first the TPT expressions for 
these quantities. These expressions involve the backward 
committor function, q-{z), which gives the probability 
that a trajectory observed at point z is coming from A 
last rather than from B (i.e. that it is assigned to A 
rather than B using the jargon introduced above). By 
definition q-{z) = \ li z £ A, q^^z) = Q ii z ^ B, 
and < q-{z) < 1 otherwise. The backward committor 
function is useful because by the Markovian assumption 
it follows that, at statistical steady state, the probability 
density to observe a trajectory at point 2; g at any 
given time t and that this trajectory is assigned to A 
at that time is simply qa{z) = q{z) q-{z). Similarly the 
probability density to observe a trajectory at point 2; S f2 
at time t and that this trajectory is assigned to B at 
that time is qb{z) — l~ QAi^) — g{z) (1 — q-{z)). Since 
PA = QAiz)dz and pB — Jq^ QBiz)dz by definition, 
this implies that 



PA^ giz)q_{z)dz <l, Pb^I 
Jn 



PA- 



(12) 



Similarly, it is easy to see that vr is the total probability 
flux associated with gA{z) = g{z)q-{z) going through 
any dividing surface between A and B (like e.g. the 
boundary oi B). The explicit form of this flux depends 
on the specifics of the dynamics and it is given in Ref. 0: 
let us omit to repeat this formula here since it will not 
be important in what follows. 




5, = fi 



FIG. 3: Schematic illustration of a piece of trajectory crossing 
the cells Ba with a = 1, . . . , 5, where cell Bi is identified with 
the reactant state A and cell B5 with the product state B. 
The circles represent the entry points which would be used 
to perform sampling restricted in the cells. If only the red 
points are used as re-entry points, then the restricted sam- 
pling selects the pieces of this trajectory when it is assigned 
to A (these pieces are depicted in red). 

In practice, we do not know explicitly g_ (z) (nor even 
g{z) in most nonequilibrium systems) but we can still 
make use of the observations above to modify the sam- 
pling procedure explained before. Suppose that we de- 
fine the reactant state A as being the union of a group 
of cells Ba and the product state B as the union of an- 
other group. Clearly, what we would then like to do is 
modify the sampling procedure in such a way that only 
the re-entry points associated with the trajectory when 
it is assigned to A are put in the databanks (see the 
illustration in Fig. [3]) and keep track of the associated 
probability fluxes through the boundary of the cells. In- 
deed, using these re-entry points and these fluxes only, 
we would then simulate in the cells pieces of trajectories 
that are statistically indistinguishable from the unbiased 
trajectory when it is assigned to A. We could then com- 
pute the total probability in the cells to get an estimate of 
Pa (and hence ps = 1 — pa) as well as the total flux into 
B to get i/ji. We claim that there is a simple procedure 
to do these operations in practice. The key observation 
is that trajectories assigned to A can be propagated like 
regular trajectories while they are outside of A and B - 
the only constraint imposed on them is via a boundary 
condition in their past: they need to have come from A 
rather than B last. But we know what this boundary 
condition entails, at least in a statistical sense: indeed 
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the probability density gA{z) that a trajectory be at z 
and be assigned to A is equal to the statistical steady 
state probability density g{z) ioi z ^ A while it is equal 
to for z G B. We also know that the probability flux 
out of A is the statistical steady state one and the prob- 
ability flux out of B is identically zero. 

We can easily impose these boundary conditions in 
practice by modifying our sampling procedure as follows. 
Suppose that we have performed a sampling as before 
and computed the steady state (unbiased) tTq and i'a,(3- 
We can then run another independent sampling where 
we only consider the cells outside of A and B. In these 
cells, we run trajectories as before (though, as we will 
see, using re-entry points that are different from the ones 
calculated before), store their exit points from the cells 
to build databanks of re-entry points in other cells, and 
compute 



TT^ = Jg gA{z)dz which, from (fT2|) . means that 



(13) 



where N^^ is the total number of times the trajectories 
hit the boundary between B^ and Bfj and is the total 
simulation time in cell B^. Note that z^^^ is only defined 
by ([13]) if a is the index of a cell Ba not inside A and B 
(the index (3, on the other hand, runs over all the cells, 
including those forming A and B). Consistent with the 
bias we need to impose to focus on trajectories assigned 
to A, we supplement this by i'^^ — Va,f3 if ct is the index 
of a cell Ba used to define A (since the effective rate 
of exit out of A must be the unbiased statistical steady 
one) and by z^^^ = if a is the index of a cell Ba used 
to define B (since the effective rate of exit out of B must 
be zero) . We also set tt-^ = tTq, if a is the index of a cell 
Ba used to define A, tt^ = if a is the index of a cell Ba 
used to define B, and in all the other cells we compute 



A A 

E^A,,A \ " ^A,,A 

0=1 /3=1 



(14) 



where the index a runs over all the cells outside of A and 
B and we use as boundary conditions the values for i''^ a 
and TT^ set before when the index P is that of a cell used 
to define A or B. Finally, we compute the probability of 
re-entry on the edges of a cell Ba not inside A and B as 



dBpndB^ 



E 



iP + a). 



(15) 



I3',a 



This procedure automatically guarantees that we focus 
on the trajectories assigned to A and makes all the quan- 
tities above - the set of trajectories, the databanks, 
'^a p ^^"^ ""q ~ different from their unbiased statistical 
steady state counterparts. By construction we then have 



PA 



X] < 1, PB = l- PA- (16) 



Similarly we can compute i>j^ from the total flux into B: 



E 



a '^a,p- 



(17) 



a such that P such that 
Ba not in B Bf, in B 



Finally, we can get kA,B and kB,A from (fTO|) and be done. 

Comparison with TIS, FFS, and Markovian mileston- 
ing with Voronoi tessellations. Before illustrating our 
procedure to compute the reaction rate on a simple ex- 
ample, let us compare it to existing methods. First, 
let us point out that it is very similar in spirit to both 
rpjgii^j^ and FFSiiii^. Like TIS and FFS, our method 
is based on selecting only those trajectories which come 
from the reactant state A. The difference is in the way 
this selection is achieved. In particular, unlike in TIS 
and FFS, we do not require that the interfaces be or- 
dered monotonously, i.e. we do not need that a trajec- 
tory coming from A crosses all the preceding interfaces 
before reaching the next one. This offers more flexibility 
in the way the interfaces can be chosen. Here we did so 
using the edges of cells in a Voronoi tessellation because 
it is convenient, but the formalism above is clearly inde- 
pendent of that choice and can be applied to any type of 
interfaces. 

Regarding the relation with the Markovian mileston- 
ing procedure with Voronoi tessellation^, the main ad- 
vantage of the method proposed in this note to compute 
the rate is that it is exact and hence avoids completely 
the assumptions underlying milestoningiS,. On the other 
hand, the new procedure is also more costly since it re- 
quires not only to build databanks of re-entry points but 
also to do the sampling twice - once to get the unbi- 
ased statistical steady states quantities, and once more 
to get the reaction rate by twisting the dynamics. For 
systems at equilibrium, one may therefore prefer to use 
the Markovian milestoning procedure with Voronoi tes- 
sellation- which is cheaper. For nonequilibrium systems, 
this procedure is inapplicable (since it relies on the time- 
reversibility of the dynamics), but as a compromise, one 
could use the nonequilibrium sampling strategy for the 
unbiased system to compute the relevant quantities in 
Markovian milestoning and thereby avoid to make the 
second twisted sampling to compute the rate. Indeed, 
the formalism developed in Ref. 6 to approximate the 
dynamics by a continuous-time Markov chain does not 
rely on the dynamics being at equilibrium. For com- 
pleteness let us briefly recall the main objects in Marko- 
vian milestoning and indicate how to compute them in 
the present context. If, following Ref. @, we define the 
milestones as the common boundaries between any two 
adjacent Voronoi cells and denote these milestones by Si 
with i = 1, 2, . . . , A^, the key quantity to approximate 
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FIG. 4: Contourplot of the Mueller potential with the 36 
points (shown as grey dots) generating the Voronoi tessella- 
tion shown as grey lines. The reactant A and product B states 
are identified with cell Bi and Bse, respectively. 



the transitions between the milestones by a continuous- 
time Markov chain is the rate matrix whose off-diagonal 
elements can be estimated as 




if i?, 7^ 
if R, = 0. 



where 



a=l 



A 

a=l 



(18) 



(19) 



Here Tq, is the total simulation time in cell Ba (as before) , 
whereas A^^" is the total number of times the trajectory 
went from milestone Si to milestone Sj and Rf is the to- 
tal amount of time the trajectory is assigned to milestone 
Si in cell Ba, i.e. the total amount of time this trajectory 
is such that Si was the edge of cell Ba it hit last. We 
refer the reader to Ref. for more details. 

Illustrative example. We end this note by illustrat- 
ing our sampling procedure on the example of a system 
evolving by Langevin dynamics on a two-dimensional po- 
tential, i.e. z[t) = {x{t),v{t)), x,v eM? and 



x[t) = ■u(t), 

v{t) = -VV{x{t)) - jv{t) + y^2iF^r]{t). 



(20) 



Here V{x) is the Mueller potential^^ whose contourplot is 
shown in Fig.Hl (3 = l/{kBT) is the inverse temperature 
and r](t) is a Gaussian white-noise with mean zero and 
covariance {ili{t)Vj{t')) — SijS{t — t'). 7 is the friction co- 
efficient and, for simplicity, we have set the mass tensor to 
the identity. Below we took 7 = 100 and = 20 (which 
is about 20% the value of the energy barrier between the 



minimum at the top left corner of the Mueller poten- 
tial and the saddle point at (-0.8,0.6)). This example is 
obviously very simple and serves no other purpose than 
being a benchmark: the application of our sampling pro- 
cedure to more interesting examples will be reported else- 
where. To avoid confusions, before presenting our results 
for (|20p let us note that this system is an equilibrium one, 
with equilibrium density g{x,v) = exp{—PH{x,v)) 



where H{x, v) 



\v\ + V{x) is the Hamiltonian and 



Z = /i[j2xR2 eKp{—PH(x,v))dxdv the partition function. 
We are, however, primarily interested in computing the 
reaction rates between the reactant state A and the prod- 
uct state B shown in Fig. SI which we will achieve by 
twisting the dynamics as explained before to focus on 
trajectories assigned to A. In so doing, we automati- 
cally put the system out of equilibrium since A becomes 
a source and B a sink (and gA{x,v) ^ g{x,v)). This 
is why we need the nonequilibrium formalism developed 
above to compute reaction rates even in the case of an 
equilibrium system such as (|^n|) . 

Shown in Fig. |4] are the 36 cells that we used in our 
calculations. These cells were defined as 



Ba = {ix,v) e 



< \x 



for all f3 ^ a}, 



(21) 



where | • | denotes the Euclidean norm and Xa with 
a = 1, ... ,36 are points chosen randomly in the region 
where the energy is below a certain threshold value. We 
identified the reactant and product states with two single 
cells in the neighborhood of the two deep minima in the 
potential landscape: A = Bi and B — B^q (see FigH]). 
Note that we intentionally took many cells and disposed 
them in a way that may not be optimal for the reaction to 
check that our procedure is robust against such choices. 
Note also that by defining the cells primarily in position 
space as we do in (j2ip we guarantee that we can count 
successive crossing of their boundaries. This would not 
be the case if the boundaries of the cells were bent in ve- 
locity space, because of the noise term in ^U\\ that acts 
on the velocities. 

Both steps of the sampling procedure (the one to com- 
pute the unbiased equilibrium quantities and the other 
with the twisted dynamics) were performed as described 
above by running simulations for 10® steps in each cell us- 
ing the second order integrator of Ref. |l8| with a timestep 
At = 10^^. In the first step of the procedure we took 
TTa equal in each cell as initial condition; in the second 
step, we took — TTa as initial condition. We compared 
the results of our procedure with those obtained by gen- 
erating a long (10^° steps) unbiased trajectory by brute 
force simulation and computing vji, kA.B and ks^A from 
finite T approximations of the limits in ([5]) and 

The main outputs of our procedure are the rate of reac- 
tive trajectories uji and the reaction rates fc^.s and ks.A 
which are reported in Table |T1 These are within statis- 
tical errors of the corresponding quantities estimated by 
brute force simulation. Our procedure also produces the 
equilibrium probabilities tTq: as shown in Fig. O these 
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kA,B 


kB,A 


Twisted dynamics 
Direct simulation 


5.910"^ 
5.8 lO'^* 


7.4 10"^ 
7.1 10"^ 


3.1 10"^ 
3.2 10~^ 



TABLE I; Rate of reactive trajectories and reaction rates be- 
tween the sets A and B in the Mueller potential (see Fig. |4]) 
obtained by the procedure presented here and by direct cal- 
culation using a long unbiased trajectory. 



are within statistical errors of the corresponding values 
obtained by the brute force simulation. Also shown in 
Fig. [5] are the probabilities computed by twisting the 
dynamics. As expected, the closer to B, the smaller tt^ 
is compared to tTq. 




FIG. 5: Probability ttq to be in the cells (semi-log scale) ob- 
tained via sampling restricted in the Voronoi cells (black dots) 
and by a brute force unbiased trajectory (triangles). The 
agreement is excellent. Also shown as circles are the steady 
state nonequilibrium probabilities vr^ of the twisted dynam- 
ics. Note that by construction tv^q = (since -Bae = B) and 
we did not plot this point. 
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